Gluon saturation and the Froissart bound: a simple approach 
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At very high energies we expect that the hadronic cross sections satisfy the Froissart bound, which 
is a well-established property of the strong interactions. In this energy regime we also expect the 
formation of the Color Glass Condensate, characterized by gluon saturation and a typical momentum 
scale: the saturation scale Q a . In this paper we show that if a saturation window exists between 
the nonperturbative and perturbative regimes of Quantum Chromodynamics (QCD) , the total cross 
sections satisfy the Froissart bound. Furthermore, we show that our approach allows us to describe 
the high energy experimental data on pp/pp total cross sections. 

PACS numbers: 12.38.-t,12.38.Aw,12.38.Bx 



I. INTRODUCTION 

Understanding the behavior of high energy hadron re- 
actions from a fundamental perspective within Quantum 
Chromodynamics (QCD) is an important goal of particle 
physics. Since it was observed that the total hadronic 
cross sections grow with the center of mass energy (\/s), 
much theoretical effort has been devoted to explain this 
growth. In particular, a QCD based explanation for this 
rising behavior was proposed by Gaisser and Halzen [lj]: 
the cross section would grow because partons would start 
to play a role in the hadronic reactions. At higher ener- 
gies, lower values of the Bjorken x are accessible and the 
parton distributions (especially the gluon distribution) 
grow very rapidly leading to rising cross sections. The 
basic idea of this approach, called minijet model, is that 
the total cross section can be decomposed as follows (see 
Fig. HI (a)): 



O~tot = 




= (T + o>qod(1) 



where Co characterizes the nonperturbative contribution, 
which is in general taken as energy-independent at high 
energies 0, and <j p qcd is calculable in perturbative 
QCD with the use of an arbitrary cutoff at low trans- 
verse momenta pq. Unfortunately, this approach implies 
a power-like energy behavior for the total cross section, 
violating the Froissart bound, which is a consequence of 
the unitarity of the S matrix and states that total cross 
sections cannot grow faster than In 2 s as s — > oo. This 
bound is a well-established property of the strong inter- 
actions and puts a strict limit on the rate of growth with 
energy of the total cross sections. Over the years several 
solutions have been proposed to cure the too fast growth 
found in the minijet model most of them using the 
cikonal formalism [31. 



In parallel with these developments the study of 
the high energy limit of the linear evolution equations 
(DGLAP and BFKL) @] revealed that they should be 
modified and gluon recombination effects (consequence 
of the high density of gluons) should be included in the 
QCD evolution (5J. This expectation can be easily un- 
derstood: while for large momentum transfer k±, the 
BFKL equation predicts that the mechanism g — > gg 
populates the transverse space with a large number of 
small size gluons per unit of rapidity (the transverse size 
of a gluon with momentum k± is proportional to l/k±), 
for small k± the produced gluons overlap and fusion pro- 
cesses, gg — > g, are equally important. Currently, one 
believes that the small- a; gluons in a hadron wave func- 
tion should form a Color Glass Condensate (CGC) which 
is described by an infinite hierarchy of coupled evolution 
equations for the correlators of Wilson lines [1, 0| ■ This 
new state of matter is characterized by gluon saturation 
and by a typical momentum scale, the saturation scale 
Q s , which grows with the energy and determines the crit- 
ical line separating the linear and saturation regimes of 
the QCD dynamics. The saturation effects are small for 
k± > Qs an d very strong for k\ < Q 2 . Experimen- 
tally, there are strong evidences of nonlinear (satura- 
tion) effects at DESY-HERA. In particular, the DESY 
ep HERA data in the small- a; and low-Q 2 region can 
be successfully described in terms of saturation models 
[1, [i| [l(| HH, EH EE EH , with the measured cross sections 
presenting the geometric scaling property [l5l |. which is 
an intrinsic property of the CGC physics. Moreover, the 
CGC physics is able to describe quite well the dAu RHIC 
data (see, e.g. Ref. 16]). These results give strong sup- 
port to the existence of a saturation regime in the QCD 
dynamics (for recent reviews see, e.g., Ref. [6j). 

Some attempts to reconcile the QCD parton picture 
with the Froissart limit using saturation physics were 
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FIG. 1: Schematic behavior of the momentum distribution in 
the minijet model (a) and in the model proposed in this paper 
(b). While in the minijet model the region Aqcd < Pt < Po 
is disregarded, the region Aqcd < Pt < Qs(s) is included 
in our model and its contribution to the total cross section 
increases when the energy rises. 



firoposed in recent years, but the question remains open 
171 Il8l . Il9| | . In this paper we propose a very simple phe- 
nomenological approach to treat this problem. In the 
next section we briefly describe the minijet model and 
how we include saturation effects in it. In the subse- 
quent sections we present our numerical results and dis- 
cuss them. 



II. THE MINIJET MODEL WITH SATURATION 

In what follows we generalize the minijet model as- 
suming the existence of the saturation regime in the high 
energy limit. More precisely, we assume the existence of 
a saturation window between the nonperturbative and 
perturbative regimes of QCD, which grows when the en- 
ergy increases (since Q s grows with the energy). This 
window is shown in Fig. [1] (b). We now generalize Eq. 
((T|) introducing the saturation window: 



Otot 



2 da 

dp T — 



CO + <7sat + VpQCD , 



QCD 



A 2 da 



s/4 j 
A 2 d(T 

p diptp 
(2) 



where the saturated component, a sa t, contains the dy- 
namics of the interactions at scales lower than the satu- 
ration scale. In this region nonlinear effects are strong, 



changing the pr behavior of the differential cross section, 
which becomes much less singular in the low p? region, 
as we can see in Fig. [T] (b). The simple QCD collinear 
factorization formulas do not in general apply in this re- 
gion. 



A. CTpQCD 

At high transverse momenta we keep using the same 
collinear factorization formula employed in [l|, |3[ with 
the necessary update of the parton densities. The satu- 
ration scale arises as a natural cut-off at low transverse 
momenta of the perturbative (minijet) cross section com- 
ponent, (JpQCD, which is now given by: 



VpQCD 



I dpip / dxi dx 2 fi (xi, p%>) fj (x 2 ,Pr)^i 



(3) 



where fi{x,Q 2 ) is the parton density of the species i in 
the proton extracted from deep inelastic scattering (DIS) 
and &ij is the leading order elementary parton-parton 
cross section. At very high energies the cross section 
([3|) is dominated by gluon-gluon interactions. In what 
follows we use the MRST leading order parton distribu- 
tions [20] in our calculations of <j p qcd- Similar results 
are obtained using for instance the CTEQ6-LO parton 
distribution sets [211 ]. 

In order to evaluate <j p qcd we need to specify Q s , 
which is determined by the solution of the nonlinear evo- 
lution equation associated to CGC physics @, Q. It is 
given by: 



(4) 



where x is the Bjorken variable, with Q\ — 0.3 GeV 2 
and Xq — 0.3 x 10~ 4 fixed by the initial condition. The 
saturation exponent A has been estimated considering 
different approximations for the QCD dynamics, being 
w 0.3 at NLO accuracy [13], in agreement with the HERA 
phenomenology, where the parameters Qoj and A were 
fixed by fitting the ep HERA data [f| . 

When we go from deep inelastic scattering to hadron 
- hadron collisions, there is some ambiguity in the defi- 
nition of the equivalent of the Bjorken x. Following [23| 
we take x to be 
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(•5) 



where go is a momentum scale to be determined. From 
(HI) and (0 we immediately see that: 



Q 2 s (s)<xs> 



(G) 



with the constant of proportionality being determined by 
HERA data and our choice of qo. 
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As it will be seen, in the high energy limit <j p qcd is 
the most important contribution to a to t- With a constant 
infra-red cut-off, as in (fT|), it would grow too fast. The 
introduction of a cut-off increasing with energy can tame 
this growth, since the bulk of the integral in ([3]) comes 
from the low momentum region. This procedure was al- 
ready employed, for example in 24], in a purely prag- 
matic approach. Here we establish a connection between 
this cut-off and the energy behavior of Q s , controlled by 
the parameter A in ©. Our procedure is more physical 
and, at the same time, imposes restrictions on A. 



B. 



In order to calculate the total cross section we also 
need to specify the saturated component. There are a 
few models for o sat [H, [2^, [2|| . Most of them are formu- 
lated in the color dipole picture, in which the projectile 
proton is treated as a color dipole, which interacts with 
the target proton. We shall use the model proposed in 
Ref. [23] , in which the total cross section is given by: 



&sat 



d 2 r\%(r)\ 2 a dip {x,r) 



(J) 



where r is the dipole transverse radius and the proton 



wave function ^> p is chosen to be: 



r 2 \ 
2Sl) 



(8) 



with S p = 0.74 fm and the dipole-proton cross section 
reads: 

(T d ip(x, r ) = 2 f d 2 bAf(x, r, b) = aj\f(x, r) (9) 

where a = 2irR 2 , with R p = 0.9 fm. The dipole scatter- 
ing amplitude, Af(x,r 1 b) 1 should be given by the impact 
parameter dependent solution of a non-linear evolution 
equation, such as the Balitsky-Kovchegov equation 
A complete solution is not yet available and we would 
have to use models for Af(x, r, b). In 23] it was assumed 
that Af(x, r, b) falls exponentially with b. In this case 
Ogat does not violate the Froissart bound. However in 
[18l ] it has been argued that the dipole amplitude decays 
only as a power of the impact parameter in the periphery 
of the proton and this dependence will, after integration 
in 6, lead to logarithmic divergences. Here we prefer to 
avoid the use of models and, instead, assume the factor- 
ization implied by the second equality in @. With this 
assumption we decouple the impact parameter and en- 
ergy dependences and focus only on the energy behavior 
of the dipole and hadron-hadron cross sections. 

In the literature there are many parameterizations of 
the dipole amplitude. A brief discussion of the features 
of some recent ones can be found in [l6| . In what fol- 
lows we shall use two of them, which were shown to give 



a reasonable description of both HERA and RHIC data 
Most of the parameterizations follow the Glauber- 
like formula originally introduced by Golec-Biernat and 
Wiisthoff [l3| . The differences among them are in the 
anomalous dimension, 7. In the KKT model [9j the ex- 
pression for the quark dipole-target forward scattering 
amplitude is given by 0]: 



AA(r, x) = 1 — exp 



(10) 



where Q 2 
7(y, r 2 ) is 

l{Y,r 2 ) 



f- Q 2 and the anomalous dimension 



1 / g(F,r 2 ) \ 

2 I, + e(y, r 2 ) + V2 i{Y, r 2 ) + 7C(3)^ ) ' 

(11) 

with c a free parameter ( which was fixed in to c = 4) 
and 



£(» 2 ) 



In [l/(r 2 Q 2 )] 
(X/2)(Y-Y ) 



(12) 



The authors of Q assume that the saturation scale can 

be expressed by Q 2 {Y) = A 2 A 1 / 3 The form of the 

anomalous dimension is inspired by the analytical solu- 
tions to the BFKL equation. Namely, in the limit r — ► 
with Y fixed we recover the anomalous dimension in the 
double logarithmic approximation 7 « 1 — \Jl/{2 £). In 
another limit of large Y with r fixed, Eq. (jlip reduces 
to the expression of the anomalous dimension near the 
saddle point in the leading logarithmic approximation 
T ~ 3 14 e C(3) ' Therefore Eq. (fTTj) mimics the onset 
of the geometric scaling region [ll|, [l2[ ■ In the calcu- 
lations of Ref. it is assumed that a characteristic 
value of r is r ~ 1/(2 It?) where kx is the transverse mo- 
mentum of the valence quark and 7 was approximated 
by 7(y, r 2 ) rj 7(Y, 1/(4 kj)). In the above expressions 
the parameters A = 0.6 GcV and A = 0.3 are fixed by 
DIS data [Hi]. Moreover, the authors assume Yq = 0.6. 
The initial saturation scale used in (Tl2|) is defined by 
Q 2 Q = Q 2 (Yq) with Yo being the lowest value of rapidity 
at which the low- a; quantum evolution effects are essen- 
tial. As demonstrated in Ref. [9] this parameterization 
is able to describe the dAu RHIC data when the forward 
dipole cross section is convoluted with the respective frag- 
mentation function and the parton distributions for the 
deuteron. 

In Ref. [3] another phenomenological saturation 
model has been proposed in order to describe the dAu 
RHIC data (hereafter denoted DHJ model). The basic 
modification with respect to the KKT model is the pa- 
rameterization of the anomalous dimension which is now 
given by 



7 (F,r 2 ) = 75 . + A 7 (r,r 2 ) 



(13) 
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FIG. 2: Energy behavior of the perturbative (solid lines) and 
saturated components (dashed lines) of the total cross section 
(normalized by In 2 s and in arbitrary units) for two different 
values of the exponent A. In (a) and (b) we use the KKTm and 
GKMN dipole cross sections respectively. The perturbative 
component is the same. 



where 

A 7 (V,r 2 ) 



(1-7.) 



| log 



r 2 Qi 



XY + | log 



(14) 



with Qt — Qs(Y) a typical hard scale in the process, 
A = 0.3 and d — 1.2. Moreover, 7 S = 0.63 is the anoma- 
lous dimension of the BFKL evolution with saturation 
boundary condition. Similarly to the KKT model this 
model is able to describe the dAu RHIC data. 

As already discussed in Ref. [l(J, based on the uni- 
versality of the hadronic wave function predicted by the 
CGC formalism, we might expect that the KKT and DHJ 
parameterizations would also describe the HERA data 
on proton structure functions in the kinematical region 
where the saturation effects should be present (small x 



jepr 

and low Q 2 ). However, as shown in T6J, this is not the 
case and neither KKT nor DHJ give an acceptable de- 
scription of the HERA data on F 2 . 

Following Ref. [27| we consider a modification of the 
KKT model assuming that the saturation momentum 




1000 10000 
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FIG. 3: Energy behavior of the total pp/pp cross section. 
The nonperturbative, perturbative and saturated components 
are presented separately as well as their sum, the total cross 
section. The results are for A = 0.25. Data are from [28j], [29| ]. 
[Ij, [H and from 



scale is given by (j4| , Yq = 4.6, c = 0.2 and that the 
typical scale in the computation of £(Y, r 2 ) is the photon 
virtuality. This modified model will be called KKTm. 
We also use the modified version [l6[ of the DHJ model, 
called here GKMN, in which Q T = Q = 1.0 GeV, i. e. 
that the typical scale is energy independent. 

In order to calculate a sa t, it is also necessary to specify 
the Bjorken-x variable as in our previous calculations of 
(JpQCD ■ We use the same prescription and the same value 
for q . 



III. RESULTS 

In Fig. [2] we show in arbitrary units the energy behav- 
ior of the ratios u v qcd / In 2 s (solid lines) and a sa t / In 2 s 
(dotted lines) for two choices of A. As it can be seen, all 
curves grow slower than ln 2 s. For smaller values of A, 
such as A (— 0.1) the fall of the ratio shown in Fig. 2 
would be postponed to very high energies, y/s ~ 10 6 
GeV. Although the energy at which the behavior of 
the cross section becomes "sub-Froissart" may depend 
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FIG. 4: Energy behavior of the total pp/pp cross section for 
different values of the exponent A. Data are the same as in 
Fig. 3. 



strongly on A, one conclusion seems very robust: once A 
is finite, at some energy the growth of the cross section 
will become weaker than ln 2 s. 

After the study of the main properties of the pertur- 
bative and saturated components we can calculate the 
total cross section and compare with experimental data, 
obtained at CERN [1|, at Fermilab Tevatron HU, [U 
and in cosmic ray experiments [31, 32] . The latter refer to 
proton-air cross sections and were translated to proton- 
proton cross sections in the phenomenological study of 
Refs. [HI1. 

In Fig. [3] we show the sum ao +<J sa t + o~r>pcp compared 
with the experimental data from Refs. [H, [2i| HH, H2| • 
<To was taken to be 34 mb. In the figure the upper and 
lower panels were calculated with the KKTm and GKMN 
models, respectively. The values of go were go = 0.044 
GeV and go = 0.038 GeV respectively. 

Considering that there is only one free parameter (go) 
in our approach, we obtain a good agreement with data. 
Moreover our predictions satisfy the Froissart bound. 
Probably a better agreement may be obtained if other 
quantities are treated as free parameters, as for instance 
the effective exponent A, and included in a fitting pro- 
cedure. As already mentioned, our results have a strong 
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FIG. 5: Energy behavior of the total pp/pp cross section for 
different values of the exponent A. Data are the same as in 
Fig. 3. 

dependence on this quantity, as it can be seen in Fig. 2J 
where we present our results for the total cross section for 
different values of A. It is important to emphasize that 
there is only a small range of values of A which allow us to 
describe the experimental data. If, for instance, A = 0.4 
the resulting cross section is very flat and clearly below 
the data, while if A = 0.1 the cross section grows very 
rapidly deviating strongly from the experimental data. 
The best choice for A is in the range 0.25 — 0.3, which is 
exactly the range predicted in theoretical estimates us- 
ing CGC physics and usually obtained by the saturation 
models for the DESY ep HERA data. 

In the theory of the CGC the parameter A changes with 
the energy, being a function of the variable Y — ln(l/x). 
Since our analysis is applied to a wide range of energies 
we have included the energy dependence of A as estimated 
in [22j], which can be parameterized as: 

A = 0.3 - 0.003(V-5) (15) 

In Fig. [5] we compare the cross sections obtained with a 
fixed value of A (= 0.25) and obtained with a "running" 
A, according to (fTo*)) . As it can be seen, the difference 
between them is small. 



IV. CONCLUSIONS 

In this paper we have proposed a simple model for the 
total pp/pp cross section, which is an improvement of 
the minijet model with the inclusion of a window in the 
PT-spectrum associated to the saturation physics. Our 
model implies a natural cutoff for the perturbative calcu- 
lations which modifies the energy behavior of this com- 
ponent, so that it satisfies the Froissart bound. More- 
over, including the saturated component (calculated with 
a dipole model), we obtain a satisfactory description of 
the experimental data. Our results for the total pp/pp 
cross section also satisfy the Froissart bound. Finally, 
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we find a very interesting consistency between our model 
and the saturation models used to describe the HERA 
data: similar values of A are needed to describe both set 
of experimental data. 

In other similar approaches, such as [23] the saturated 
cross section is used over the entire pr domain, or equiv- 
alently, for dipoles of all sizes. This procedure has two 
disadvantages: it requires the introduction of an (model 
dependent) impact parameter dependence of the dipole 
cross section and it does not make use of the collinear 



factorization formula and the parton densities, which, in 
the high p? region, are very well studied both theoret- 
ically and experimentally. In this sense our work is an 
improvement on [l[ and on [23| as well. 

Acknowledgements: This work was partially financed 
by the Brazilian funding agencies FAPESP, FAPERGS 
and CNPq. 



[1] T. K. Gaisser and F. Halzen, Phys. Rev. Lett. 54, 1754 
(1985). 

[2] See, for example, H. G. Dosch, E. Ferreira and 
A. Kramer, Phys. Rev. D 50, 1992 (1994); H. G. Dosch, 
F. S. Navarra, M. Nielsen and M. Rueter, Phys. Lett. B 
466, 363 (1999). 

[3] L. Durand and H. Pi, Phys. Rev. D 40, 1436 (1989); 
Nucl. Phys. Proc. Suppl. 12, 379 (1990); R. M. God- 
bole, A. Grau, G. Pancheri and Y. N. Srivastava, Phys. 
Rev. D 72, 076001 (2005); E. G. S. Luna, A. F. Martini, 
M. J. Menon, A. Mihara and A. A. Natale, Phys. Rev. D 
72, 034019 (2005). 

[4] For a review of these equations see, for example, V. 
Barone and E. Predazzi, "High Energy Particle Diffrac- 
tion", Springer Verlag, Berlin, (2002). 

[5] L. V. Gribov. E. M. Levin and M. G. Ryskin, Phys. Rep. 
100, 1 (1983). 

[6] E. Iancu and R. Venugopalan, |arXiv:hep-p h/0303204 

A. M. Stasto, Acta Phys. Polon. B 35, 3069 (2004); 

H. Weigert, Prog. Part. Nucl. Phys. 55, 461 (2005); 

J. Jalilian-Marian and Y. V. Kovchegov, Prog. Part. 

Nucl. Phys. 56, 104 (2006). 
[7] I. Balitsky, Nucl. Phys. B 463, 99 (1996); 

Y. V. Kovchegov, Phys. Rev. D 60, 034008 (1999); 

Phys. Rev. D 61, 074018 (2000). 
[8] J. Bartels, K. Golec-Biernat and H. Kowalski, Phys. Rev 

D66, 014001 (2002); H. Kowalski and D. Teaney, Phys. 

Rev. D 68, 114005 (2003). 
[9] D. Kharzeev, Y.V. Kovchegov and K. Tuchin, Phys. Lett. 

B599, 23 (2004). 
[10] M. S. Kugeratski, V. P. Goncalves and F. S. Navarra, 

Eur. Phys. J. C 44, 577 (2005). 
[11] E. Iancu, K. Itakura, S. Munier, Phys. Lett. B590, 199 

(2004). 

[12] E. Iancu, K. Itakura and L. McLerran, Nucl. Phys A708, 
327 (2002). 

[13] K. Golec-Biernat and M. Wiisthoff, Phys. Rev. D 59 

(1999) 014017, ibid. D60, 114023 (1999). 
[14] A. Dumitru, A. Hayashigaki and J. Jalilian-Marian, Nucl. 

Phys. A 765, 464 (2006); Nucl. Phys. A 770, 57 (2006). 
[15] A. M. Stasto, K. Golec-Biernat and J. Kwieciiiski, 

Phys. Rev. Lett. 86, 596 (2001); V. P. Goncalves and 

M. V. T. Machado, Phys. Rev. Lett. 91, 202002 (2003); 

JHEP 0704, 028 (2007); C. Marquet and L. Schoeffel, 

Phys. Lett. B 639, 471 (2006). 
[16] V. P. Goncalves, M. S. Kugeratski, M. V. T. Machado 

and F. S. Navarra, Phys. Lett. B 643, 273 (2006). 



[17] E. Ferreiro, E. Iancu, K. Itakura and L. McLerran, Nucl. 

Phys. A 710, 373 (2002). 
[18] A. Kovner and U. A. Wiedemann, Phys. Rev. D 66, 

051502 (2002); Phys. Rev. D 66, 034031 (2002); Phys. 

Lett. B 551, 311 (2003). 
[19] T. Ikeda and L. McLerran, Nucl. Phys. A 756, 385 

(2005). 

[20] A. D. Martin, R. G. Roberts, W. J. Stirling and 
R. S. Thorne, Eur. Phys. J. C 23, 73 (2002). 

[21] J. Pumplin, D. R. Stump, J. Huston, H. L. Lai, P. Nadol- 
sky and W. K. Tung, JHEP 0207, 012 (2002). 

[22] D. N. Triantafyllopoulos, Nucl. Phys. B 648, 293 (2003). 

[23] J. Bartels, E. Gotsman, E. Levin, M. Lublinsky and 
U. Maor, Phys. Lett. B 556, 114 (2003). 

[24] K. J. Eskola, K. Kajantie and K. Tuominen, Nucl. Phys. 
A 700, 509 (2002); K. J. Eskola, P. V. Ruuskanen, 
S. S. Rasanen and K. Tuominen, Nucl. Phys. A 696, 
715 (2001). 

[25] A. I. Shoshi, F. D. Steffen and H. J. Pirner, Nucl. Phys. 

A 709, 131 (2002). 
[26] B. Z. Kopeliovich, I. K. Potashnikova, B. Povh and 

E. Predazzi, Phys. Rev. Lett. 85, 507 (2000); Phys. Rev. 

D 63, 054001 (2001). 
[27] M. V. T. Machado, Eur. Phys. J. C 47, 365 (2006). 
[28] A. S. Carrol et al., Phys. Lett. B61, 303 (1976); A. S. 

Carrol et al., Phys. Lett. B80, 423 (1979); U. Amaldi, 

and K. R. Schubert, Nucl. Phys. B166, 301 (1980); V. 

Bartenev et al., Phys. Rev. Lett. 31, 1089 (1973); N. A. 

Amos et al., Nucl. Phys. B262, 689 (1985); A. Breakstone 

et al., Nucl. Phys. B248, 253 (1984); M. Ambrosio et al., 

Phys. Lett. B115, 495 (1982); C. Augier et al., Phys. Lett. 

B316, 448 (1993). 
[29] C. Augier et al., Phys. Lett. B 344, 451 (1995); M. Bozzo 

et al., Phys. Lett. B 147, 392 (1984); F. Abe et al., Phys. 

Rev. D 50, 5550 (1994); N. A. Amos et al., Phys. Rev. 

Lett. 68, 2433 (1992); C. Avila et al., Phys. Lett. B 537, 

41 (2002). 

[30] C. Albajar et al., Nucl. Phys. B 309, 405 (1988). 

[31] Fly's Eye Collab., R.M. Baltrusaitis et al., Phys. Rev. 

Lett. 52, 1380 (1984). 
[32] AGASA Collab., M. Honda et al., Phys. Rev. Lett. 70, 

525 (1993). 

[33] M.M. Block, F. Halzen and T. Stanev, Phys. Rev. D 62, 

077501 (2000). 
[34] N. N. Nikolaev, Phys. Rev. D 48, 1904 (1993). 



